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In order to verify Percival's conjecture [J. Phys. B 6,L229 (1973)] we study a planar billiard in 
its classical and quantum versions. We provide an evaluation of the nearest-neighbor level-spacing 
distribution for the Cassini oval billiard, taking into account relations with classical results. The 
statistical behavior of integrable and ergodic systems has been extensively confirmed numerically, but 
that is not the case for the transition between these two extremes. Our system's classical dynamics 
undergoes a transition from integrability to chaos by varying a shape parameter. This feature allows 
us to investigate the spectral fluctuations, comparing numerical results with semiclassical predictions 
founded on Percival's conjecture. We obtain good global agreement with those predictions, in 
clear contrast with similar comparisons for other systems found in the literature. The structure of 
some eigenfunctions, displayed in the quantum Poincare section, provides a clear explanation of the 
conjecture. 
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I. INTRODUCTION 



In 1973 Percival conjectured that in the semiclassical limit, the spectrum of a generic dynamical system consists of 
two parts with strongly contrasting properties: a regular and an irregular part 0. At the classical level such a system 
exhibits a mixed dynamics: regular regions dominated by tori and chaotic regions with mixing behavior coexist in 
the phase space. 

In order to characterize a semiclassical spectrum it is advantageous to consider a sensitive fluctuation measure. 
The probability distribution p(s) of the spacing s between succesive levels is of particular interest because it contains 
information of the spectrum on its finest scale. In the special case of multidimensional integrable systems, Berry and 
Tabor || showed that the levels are uncorrelated and p(s) is governed by a Poisson distribution. The other special 
case corresponds to mixing systems where almost all orbits explore densely and chaoticaly the energy surface. In this 
case, Bohigas et al. || conjectured that the fluctuation properties of these spectra can be modeled by the ensemble of 
random real symmetric matrices [the Gaussian Orthogonal Ensemble (GOE)] For such systems, where the energy 
levels display repulsion, p(s) is closely approximated by the Wigner distribution; however, we have used the exact 
distribution ^ because the differences are meaningful, as we shall see below. 

In the generic case, Berry and Robnik ||, based on Percival's conjecture, considered independent sequences of 
levels associated with each connected regular or irregular classical phase-space region. When only one chaotic region 
predominates (the situation considered in this article), they provided an expression for the distribution p(s) in terms 
of the classical fraction p cl of regular regions [referred to as the Berry- Robnik distribution (BRD)]. 

Although the special cases have been extensively confirmed numerically [j7l [9]J^, [lo| , |lll , with some well understood 
exceptions [ ^2|Jl4| , for generic systems numerical calculations give rise to contradictory conclusions. Recently, there 
has been a number of numerical works |]l5|--|l7| showing that the Brody distribution (BD) jL8) gives quite a satisfactory 
fit globally. The Brody distribution is a one-parameter family of distributions that interpolates between Poisson and 
Wigner in a simple way; however, it has no semiclassical meaning. On the other hand, in 1994 Prosen and Robnik 
confirmed numerically semiclassical predictions (the BRD) working on an abstract dynamical system : the standard 
map on a torus. To agree with this theory, they needed to compute extremely high excited states (around the 30 
10 6 ). However, at not so excited states they found good global agreement with the Brody distribution both for the 
standard map and for the Limacon-Robnik billiard. In 1995 Prosen arrived at the same conclusion working on a 
two-dimensional semiseparable oscillator pCj ], 

The Brody-like behavior at small spacing is understood in terms of tunneling between classically separated regions 
of phase space; however, the global agreement with the Brody distribution has no theoretical foundations. On the 
other hand, the very slow trasition to semiclassical predictions can be explained by the presence of partial barriers in 
the chaotic regions because the corresponding statistics is not a GOE for finite Ti Jl(| 13 . 

The goal of the present article is to verify that the classical support of eigenfunctions can be clearly identified as 
regular or chaotic and this classification is only affected by tunneling between classically separated regions of phase 
space (this effect decays exponentially when h decreases, as it was pointed out in ref. ||). We compute the spectral 
fluctuations of a one-parameter family of planar billiards : the Cassini ovals |2If| . The classical dynamics in this 
billiard is mixed, going from integrability to chaos by varying a single parameter. We have chosen this parameter 
such that the classical dynamics does not show partial barriers inmerse in the chaotic sea to study the Percival's 
conjecture through the accuracy of the BRD. Moreover, we study qualitatively the eigenfunctions in phase space to 
provide additional support to the results. 

Our work is organized in the following way. In Sec. || we introduce the classical system. Section [II is devoted 



to the description of the quantum system (its energy spectrum and the corresponding eigenfunctions). In Sec. IV we 
study the resulting energy level statistic. Finally, Sec. [V]is devoted to conclusions. 



II. THE CASSINI OVAL BILLIARD 



Our billiard consists of a free-moving point particle inside a two-dimensional box that bounces off the boundary 
elastically. The boundary of our billiard system is given by a fourth-order curve, the Cassini oval: 

rir 2 = a 2 , 

where r\ and r 2 are distances from two foci located at x = ±c and y — 0. In cartesian coordinates it can be cast into 
the form: 

(x 2 + y 2 ) 2 - 2c 2 (x 2 ~ y 2 ) = a 4 - c 4 . (2.1) 
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We have two characteristic lengths. However, the shape of the boundary is defined by the ratio d — a/c (from now on, 
the shape parameter), which determines the following boundary types: y/2 < d, the boundary is an oval; 1 < d < y/2, 
the boundary is an oval with a neck; d < 1, the boundary becomes unlinked (two ovals separate). 

In the present work we investigate d > 1 values. Decreasing the shape parameter, the classical behavior goes from 
the regular motion (when d — > oo, the boundary is a circle) to the chaotic one (when d — > 1). Using the reflection 
symmetries of the boundary, we consider the motion in the region x > 0,y > (desymmetrized billiard). That is, we 



study the quarter billiard defined by the boundary (eq. 2.1) for x > 0,y > and the coordinate axes x and y. We 



study the Cassini oval billiard for two values of d. One of them is d = \/2 (the value of the parameter for which the 
neck begins to appear; see Fig. |l|(a)) whose form mimics the Bunimovich stadium billiard p2[. Figure [2|( a) shows the 
Poincare surface at the boundary using Birkhoff coordinates. The coordinate q is related to the arclength coordinate 
at the boundary where the bounce takes place by q = (arclength)/ (perimeter); and p = p • t/|p| is the fraction of 
tangential momentum at this point (t being the tangent unit vector to the boundary). Exploiting the time-reversal 
symmetry we show only the p > 0, < q < 1/4 region. The classical phase space has a bouncing-ball regular region 
dominated by invariant curves and a chaotic region with unstable short periodic orbits equivalent to those appearing 
in the Bunimovich stadium billiard. A resonance of winding number 6 defines the last great regular region before 
chaos begin to appear (an eigenfunction existing on the chain of islands defined by the resonance is shown in Fig. 
H (a)). We have found two very small stable regions corresponding to a stable bifurcation of the unstable bow tie 
periodic orbit (two small dark dots can be observed in the chaotic region). (Figure |^ (b) shows an eigenfunction 
existing in that region of the phase space). We have not detected any other regular region embedded in the chaotic 
sea. 

The other shape we have studied (d = 2) is closer to an ellipse (see Fig. 0(b)). In this case the bouncing-ball region 
of the phase space is greater than before as it can be seen in Fig. |^(b). Moreover, a regular region appears as a thin 

band for p > 0.9 values, dominated by whispering gallery trajectories. Phase space is very mixed and many stable 
islands of very different sizes are interspersed with the chaotic trajectories. 

By selecting two regions corresponding to chaotic motion in the phase space, we have calculated diffusion times 
between them. The results for d = \/2 are independent of the chosen regions. In the other case (d = 2) this time 
is strongly dependent on them, and we have obtained diffussion times one order of magnitude greater than those of 
d = \/2. This fact is related to partial barriers such as those shown in Fig. ||(a) . 

We have determined p cl , the fraction of the phase space that corresponds to regular motion for both the values of 
the parameter. We have found that p cl = 0.172 when d = \/2 while p cl = 0.394 for d — 2. 

III. THE QUANTUM BILLIARD 

To study the quantum billiard we solve the timc-indcpcndcnt Schrodinger equation for one particle inside a two- 
dimensional box T> with Dirichlct boundary conditions at the impenetrable walls dT>: 

V 2 0(r) = ~k 2 <j){f), in V, 
<j> = 0, on dV , 

where 



k = ; 



and V corresponds to the surface of the desymmetrized billiard, that is the surface bounded by the Cassini oval (2.1) 
and the x > and y > semiaxis. So, the Dirichlet boundary condition implies that only odd-odd solutions of 
the full billiard will be found. We have employed a new technique, the scaling method p3| : this is a very efficient 
one dimensional method developed to compute eigenvalues and eigenfunctions of quite general planar billiards (for 
three dimensional billiards this is practically the only available method to obtain high excited states [p4j ). The great 
advantage of the method is that all eigenvalues and eigenfunctions in a narrow k interval are computed simultaneously 
with comparable accuracy, thus avoiding time consuming searches and the posibility of missing some state. 

We have calculated the energy levels from the fundamental state up to the 25000th level for d = y/2 and up to the 
10000th for d = 2. Moreover, we took a sequence of 5000 levels between the 62210th and 67210th for both values of 
d. We obtained the eigenvalues with an average precision of 10 -6 of the mean level spacing for chaotic eigenfunctions. 
Regular states are only limited by the computer (double) precision. 
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We have studied the eigenfunctions of the billiard in different regions of the spectrum. In general, it is possible 
to identify each eigenfunction with a classical region. This identification of the eigenfunctions is more clearly seen 
in the stellar representation |2{|. In it, the Husimi distribution of the normal derivative on the boundary represents 
the eigenfunctions in the Poincare section in Birkhoff coordinates. The stellar representation of an eigenfunction is 
compared directly with the classical Poincare section. As an example we show some eigenfunctions for d — \/2 (we 
have taken the area of the desymmetrized billiard equal to 7r/4). Fig. || (a) shows a linear density plot of the square 
of a state existing on the chain of islands defined by a resonance of winding number 6 and Fig. |] (a) shows the same 
state in stellar representation. Figs. || (b) and [| (b) show a scar of the bow tie unstable periodic orbit. Figs. || 
(c) and |] (c) show a delocalized state distributed over all the classically chaotic region (chaotic state). Finally, Figs. 
H (d) and || (d) show a bouncing-ball state. 

We stress that the states appearing in Figs. | (c) and (d) are quasi degenerate. The distance between them is 
a very small fraction (s = 0.00018) of the mean level spacing. It is clear from Figures [| (c) and (d) that these 
states practically do not interact because they exist in different regions of phase space. On the other hand, this 
example shows us that it is necessary to evaluate the eigenvalues with high accuracy in order to calculate the spectral 
fluctuations of the system (see the next section). 

For d = 2, the qualitative description of the eigenfunctions is equivalent to the previous one. However, there is a 
significant fraction of localized eigenfunctions existing in the chaotic region. One of them is shown in Fig. |^(b) and 
(c). 



IV. THE ENERGY-LEVEL STATISTIC 



In this section we analyze the level spacing distribution of the numerical data described previously. The counting 
function N(k) gives us the number of levels with wave number below k. Weyl's formula with border, angle and 
curvature corrections provides a good estimate for the smooth part < N(k) > of the counting function. 

In order to verify that no levels had been lost we have compared < N(k) > with a smoothed version of N(k). 
Defining a new sequence by K n =< N(k n ) > ( where k n belongs to the original sequence of wave numbers ) we take 
into account of the "unfolding" procedure by which a unit mean spacing is given to the series of levels llF 



Following ref. 
rather than p(s) 



in the analysis of the data, we use the cumulative level spacing distribution VT(s) = J Q p(y)dy 
because it is numerically easier to evaluate. We fit numerical curves with the analytical expression 



for the Brody family of cumulative spacing distributions: 

Wf{a) = 1 - cxp {-bs p+1 ) , 

with b = (r((/3 + 2)/ {13 + l))) /3+1 and the theoretical Berry-Robnik distribution: 

W? R (a) = 1 - cxp (-ps)[pQ GO E(s) - (1 - p){W GO e{s) - 1)] 

with Qgoe(s) = (1 — p) J s °° dX J x PGOE(y)dy, Pgoe(d) being the exact GOE spacing distribution, and Wgoe — 
1 + dQGOE/d((l — p)s). Qgoe and Wgoe are tabulated as \P and F (taking as mean density 1 — p) respectively 
in reference for instance. This exact GOE evaluation allows us to distinguish deviations of numerical data from 
theory without including the difference between Wigner and exact GOE formulas (which are of approximately the 
same order; this fact can be verified in Fig. ^| (a) which shows the difference between the BRD using Wigner surmise 
and the exact GOE results). 

Deviations of numerical data from best-fitting curves can be better seen with a transformation defined by: 



U{W) = (1/tt) cos" 1 y/(l-W) (4.1) 

which has constant statistical error over all s. Furthermore, if we plot U(W) versus W¥ we will have an equally 
spaced distribution of points on the abscissa (see Fig. ||, for example). 

We have evaluated x 2 = E^Li W s i)n«m - W{s,) theo ) 2 weighted with {8W{ Sl )) 2 = W{ Si ){\ - W{s,))/N so that 
we could find the optimal values of f3 and p BR (p BR is the resulting best fitting value for the fraction p of regular 
levels employed as a free parameter in order to find the lowest % 2 ). Results from fittings can be seen in Table |. 

It is clear that, for d — v2, BRD curve fits data much better than BD curve. Taking the sequence of 25000 first 
levels we found that \ 2 f° r the BRD was approximately five times lower than for the BD. Even for the small s range 
we could verify a better agreement between data and BRD formula than for the BD one, although deviations due to 
tunnelling can be clearly seen in Figure ^| (a). For large values of s, differences between data and BRD could evidence 
discrepancies with respect to GOE behaviour. Taking the sequence of 5000 levels between number 62210 and 67210 



4 



we got a x 2 f° r BRD that was 25 times lower than for BD. This can be checked with Fig. || (b), where the region of 
small s still shows some deviations from best fitting BRD but the range where this happens has become very narrow. 
In fact, for all values of s, the agreement is excellent and we can say that BRD is working perfectly well at these not 
so high energy levels. 

In the case of d = 2, for the first 10000 levels, a same order \ 2 was found for the two distributions although was 
better for the BRD than for the BD. In Fig. [?| (a) we can appreciate that BD fits numerical data in an acceptable 
way for small s values. This is due to tunnelling effects that persist in a wider s range than for the d — \/2 case. A 
more complicated structure of classical phase space makes tunnelling processes more significant. However, for greater 
s BD no longer follows numerical data, so there is no global agreement with it. As in the previous case, we take 5000 
levels between 62210 and 67210. For these energies, \ 2 f° r BD is twice the value for BRD. In Fig. ^ (b) we can see 
that discrepancies in the small s region have reduced, and, though not as clear as in d — \/2 case, numerical data is 
well adjusted by BRD globally. 

Finally, we investigated the behaviour of data while going from low to high energies. In order to do so we took first 
3000 levels for both of the shapes; then three stripes of levels, the first, second and third 8000 levels for d — \[2 and 
the first, second and third 4000 levels for d = 2. Results of fittings can be seen in Table [n| It is clear from Figures ^ 
and |^, that there is no transition from Brody to Berry-Robnik regime, neither for d — \/2, nor for d = 2. 

V. SUMMARY AND CONCLUSIONS 

To verify Percival's conjecture, we have studied the quantum version of a billiard, depending on one shape parameter 
d. This system shows mixed classical dynamics, going from integrability (d — > oo) to chaos (d — + 1) as the parameter 
is varied. 

On the one hand, the Husimi distribution of the eigenfuntions that we have obtained clearly displays the classical 
structure of the phase space. Though mixings among regular wavefunctions and irregular ones are expected (based on 
the abscence of degeneracies in this one-parameter system that has been desymmetrized) they seem to happen only 
when the energy difference is surprisingly small ( even for levels as low as N ~ 2500 such as the ones exemplified in 
Fig. | (c) and (d)). So, we can say that Percival's conjecture is effectively working. Mixed functions are exceptions 
mainly originated in the states whose Husimi distributions localize on the last KAM tori. 

On the other hand, the level spacing statistics was fitted by two distributions depending on one parameter: the 
semiclassical Berry-Robnik distribution (BRD), which is founded on Percival's conjecture, and the Brody distribution 
(BD). We have found, in all cases we have analyzed that the BRD is the best one. Moreover, as expected, fits of 
the BRD are better for increasing values of the wave number k; however, discrepancies due to tunnelling effects are 
found for small values of s, even for the best-fitting case (see Fig. || (b) - This subject is currently under study at the 
moment). For d = y/2, the parameters p of the BRD corresponding to best fits, do not show significant differences 
from the fraction of the regular region of the classical phase space. In the case of d = 2, the best p's, are systematically 
greater than the classical value. 

Our calculation of diffusion times shows that the cassini billiard with d = \/2 has only one chaotic region. There 
are not partial barriers dividing chaotic regions of comparable sizes. So we expect that the chaotic component of the 
statistic is given by a single GOE distribution. In the case d = 2 we have determined partial barriers between a main 
chaotic region and a little one near the whispering gallery region. The size of the latter is so little that it enlarges (for 
low energies) the regular region of the phase space rather than contributes as an independent GOE. This fact can be 
seen in Fig. |5j(b) and (c) , where we show an eigenfunction localized in this small classical chaotic region but having 
regular characteristics. Moreover, there are small regular islands immerse in the chaotic sea that are more relevant for 
d = 2 (see Fig. || ) than for d = y/2. They produce quantum localization in their chaotic neighborhood (for instance, 
see the Poincare surface section Fig. |^ (a) and the Husimi distribution Fig. ^ (b) ) and consequently the resulting p 
of the fits, overestimates the fraction corresponding to regular motion in the classical phase space. We stress that, in 
all cases, the BD has significant differences from the numerically calculated distribution. So we do not observe the 
BD to BRD transition that was seen in other systems. 
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FIG. 1. Desymmetrized Cassini oval (upper-right quarter of the curve with segments of the coordinate axes x and y). a) 
d — \pl with Birkhoff coordinates p and q. b) d = 2. 

FIG. 2. Poincare surface of section expressed in Birkhoff coordinates (p and q) for p > and < q < 0.25. a) d = \/2. b) 
d = 2. 

FIG. 3. Wave functions of Cassini oval billiard for d — \/2 with wave numbers k displayed below each one of them. We can 
see: a) regular eigenfunction, b) eigenfunction strongly localized on a scar reminiscent from an unstable periodic orbit (bow 
tie), c) irregular wavefunction and d) bouncing-ball wavefunction extremely close in energy to the previous one but without 
exhibiting mixture. 

FIG. 4. Husimi plots corresponding to the wave functions displayed in Figure 3. This shows: a) great localization on a 
classical resonance, b) remarkable localization on the classical region where an unstable orbit is found to exist, c) almost 
uniform extension over chaotic region and d) strong localization on one torus belonging to the interior of the principal regular 
island situated at the origin. 

FIG. 5. Effects of barriers: a) Poincare section for d = 2 taking only one initial condition very close to the whispering gallery 

region (p > 0.9). We can observe two partial barriers limiting flux and defining two small chaotic regions (which are filled more 
densely than the major chaotic part of phase space), b) Husimi distribution of an eigenfunction localized over one of these 
small regions, c) Same eigenfunction in configuration space. 

FIG. 6. Differences between numerical U and U BR taking the best-fitting value p BR are displayed with a solid (fluctuating) 
line for d — \/2. Dotted lines that follow the solid one represent the SU uncertainty band, a) N between 1 and 25000 and b) N 
between 62210 and 67210. In a), difference between U BR calculated with Wigner surmise and with exact GOE results is also 
plotted with a dashed line. We can see that deviations of numerical data from this curve are of the same order. 

FIG. 7. Same as in Figure 6 but for d = 2. In a) N between 1 and 10000 and b) N between 62210 and 67210 (we use the 
same line patterns as in Figure 5). 

FIG. 8. Same as in Figure 6, but for d = \/2. The levels used for plottings are: a) 1-3000, b) 1-8000, c) 8001-16000 and d) 
16001-24000. 



FIG. 9. Same as in Figure 7, but for d = 2 . Levels used in this case are: a) 1-3000, b) 1-4000, c) 3000-7000 and d) 
6000-10000. 
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TABLE I. Results of fittings. Comparison between regular fraction of levels p BR and integrable part of classical phase space 



d 


N 




X 2 /N 




X 2 /N 


P d 


V2 


1 - 25000 


0.570 


19.61 


0.177 


4.34 


0.172 


V2 


62210 - 67210 


0.599 


4.50 


0.165 


0.20 


0.172 


2 


1 - 10000 


0.230 


5.63 


0.433 


4.79 


0.394 


2 


62210 - 67210 


0.226 


2.27 


0.437 


0.98 


0.394 



TABLE II. Results of fittings. Three stripes of increasing energy are displayed for the two values of shape parameter. 



d 


N 


p b 


X"/N 


P UH 




V2 


1 - 8000 


0.59 


10.42 


0.17 


3.15 


V2 


8001 - 16000 


0.57 


5.79 


0.17 


1.33 


V2 


16001 - 24000 


0.54 


4.65 


0.19 


1.32 


2 


1 - 4000 


0.24 


3.00 


0.42 


2.88 


2 


3000 - 7000 


0.22 


2.68 


0.44 


1.02 


2 


6000 - 10000 


0.22 


1.61 


0.45 


1.97 
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